#install.packages('ggtree')
#install.packages('ggplot2')
#install.packages('dplyr')
#install.packages('treeio')
#install.packages('phytools')
#install.packages('ape')
#install.packages('ggpubr')
library(ggtree)
## Registered S3 method overwritten by 'treeio':
## method from
## root.phylo ape
## ggtree v2.2.1 For help: https://yulab-smu.github.io/treedata-book/
##
## If you use ggtree in published research, please cite the most appropriate paper(s):
##
## [36m-[39m Guangchuang Yu. Using ggtree to visualize data on tree-like structures. Current Protocols in Bioinformatics, 2020, 69:e96. doi:10.1002/cpbi.96
## [36m-[39m Guangchuang Yu, Tommy Tsan-Yuk Lam, Huachen Zhu, Yi Guan. Two methods for mapping and visualizing associated data on phylogeny using ggtree. Molecular Biology and Evolution 2018, 35(12):3041-3043. doi:10.1093/molbev/msy194
## [36m-[39m Guangchuang Yu, David Smith, Huachen Zhu, Yi Guan, Tommy Tsan-Yuk Lam. ggtree: an R package for visualization and annotation of phylogenetic trees with their covariates and other associated data. Methods in Ecology and Evolution 2017, 8(1):28-36. doi:10.1111/2041-210X.12628
library(ggplot2)
library(dplyr)
##
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
##
## filter, lag
## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, union
library(treeio)
## treeio v1.12.0 For help: https://yulab-smu.github.io/treedata-book/
##
## If you use treeio in published research, please cite:
##
## LG Wang, TTY Lam, S Xu, Z Dai, L Zhou, T Feng, P Guo, CW Dunn, BR Jones, T Bradley, H Zhu, Y Guan, Y Jiang, G Yu. treeio: an R package for phylogenetic tree input and output with richly annotated and associated data. Molecular Biology and Evolution 2020, 37(2):599-603. doi: 10.1093/molbev/msz240
library(phytools)
## Loading required package: ape
##
## Attaching package: 'ape'
## The following object is masked from 'package:treeio':
##
## drop.tip
## The following object is masked from 'package:ggtree':
##
## rotate
## Loading required package: maps
##
## Attaching package: 'phytools'
## The following object is masked from 'package:treeio':
##
## read.newick
library(ape)
library(ggpubr)
##
## Attaching package: 'ggpubr'
## The following object is masked from 'package:ape':
##
## rotate
## The following object is masked from 'package:ggtree':
##
## rotate
setwd('/Users/owlex/Dropbox/Documents/Northwestern/rcs_consult/r_phylogenetics_workshop/r_phylogenetics_worshop') # change this to your R Markdown's file path
We will learn about reading and plotting trees using a set of 20 genes that belong to four different classes of enoyl reductases (five per class) that are involved in fatty acid synthesis. These are fabK, fabG, fabV, and fabI.
For your own interpretation, know that FabV, FabK, and FabI have the same enzymatic function, but fabG does not. Despite the overlap in function, fabG and FabI are more closely related to each other by sequence. FabV is the next most closely related, and fabK is the least related.
Throughout the workshop, I highlight when I’m using non-ggtree packages by calling on the package name first and then the function. For example using the ape package’s function read, I call on it by typing ape::read(). Hopefully this makes the roles of the different packages being used more obvious and how they fit together more clear.
There are a variety of different phylogenetic tree file types. One of the most popular ones is the .newick (also abbreviated .nwk)
We can read in the file using ape’s read.newick() function.
tree_nwk <- ape::read.tree('raw_enoyl_seqs.nwk')
tree_nwk is of class ‘phylo’. When can always confirm this by using the class() function.
class(tree_nwk)
## [1] "phylo"
We can access the information about the tree object
# general info
tree_nwk
##
## Phylogenetic tree with 20 tips and 18 internal nodes.
##
## Tip labels:
## tr|B0BP34|B0BP34_ACTPJ, tr|A0A1R0ED53|A0A1R0ED53_HAEPA, tr|A0A0M3G239|A0A0M3G239_HAEHA, sp|P0AEK4|FABI_ECOLI, sp|O24990|FABI_HELPY, sp|I6Y778|FABG4_MYCTU, ...
## Node labels:
## , 0.197, 0.679, 0.990, 0.995, 0.183, ...
##
## Unrooted; includes branch lengths.
# tip label names
tree_nwk$tip.label
## [1] "tr|B0BP34|B0BP34_ACTPJ" "tr|A0A1R0ED53|A0A1R0ED53_HAEPA"
## [3] "tr|A0A0M3G239|A0A0M3G239_HAEHA" "sp|P0AEK4|FABI_ECOLI"
## [5] "sp|O24990|FABI_HELPY" "sp|I6Y778|FABG4_MYCTU"
## [7] "sp|P0AEK2|FABG_ECOLI" "sp|P9WGT3|MABA_MYCTU"
## [9] "sp|P71534|MABA_MYCS2" "sp|P24182|ACCC_ECOLI"
## [11] "sp|Q73Q47|FABV_TREDE" "sp|Q9HZP8|FABV_PSEAE"
## [13] "sp|Q62L02|FABV_BURMA" "sp|Q9KRA3|FABV1_VIBCH"
## [15] "sp|Q2P9J6|FABV_XANOM" "sp|P54616|FABI_BACSU"
## [17] "sp|Q6GI75|FABI_STAAR" "sp|Q81GI3|FABI_BACCR"
## [19] "tr|E6KYQ2|E6KYQ2_9PAST" "tr|A0A4Y9GUQ1|A0A4Y9GUQ1_9PAST"
The ggtree package makes it very easy to quickly read in a .newick file and view a nicely formatted phylogenetic tree.
tree_nwk <- ape::read.tree('raw_enoyl_seqs.nwk')
ggtree(tree_nwk)+ # read in tree
geom_tiplab()+ # show tip labels
geom_treescale(5) # recenter the tree
## Warning: `tbl_df()` is deprecated as of dplyr 1.0.0.
## Please use `tibble::as_tibble()` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_warnings()` to see where this warning was generated.
Not all phylogenetic tree files can be read with read with ape::read.tree(). The treeio package can open virtually all other existing tree file types.. Click on the provided link see their commands.
Another common tree filetype are those produced by the program RAxML. These are a modified newick file that contains bootstrap numbers at each node.
We will use the ape::read.tree() to open a raxml file in R.
Use ape::read.tree() to read in the file raw_enoyl_seqs.nwk and assign it to the variable my_nwk
Display the my_nwk tree using ggtree()
Show the geom_tiplab() for each leaf in the ggtree() visualization of my_nwk
Modify the code from part 3 so that the x-axis is rescaled at x=7 using geom_treescale()
Some file types hold additional information. For example, the RAxML tree files can contain bootstrap values at each node bifurcation.
tree_raxml <- ape::read.tree('RAxML_bipartitions.raw_enoyl_seqs')
class(tree_raxml)
## [1] "phylo"
# for a RAxML file, the bootstrap values are stored under the node.label
tree_raxml$node.label
## [1] "" "100" "54" "100" "56" "94" "47" "99" "77" "97" "100" "48"
## [13] "29" "46" "100" "85" "100" "55"
The bootstrap values can be displayed “out of the box” on the tree using the option geom_nodelabaes(label=label)
ggtree(tree_raxml)+
geom_tiplab()+
geom_nodelab(aes(label=label))+ #show bootstraps
geom_treescale(x=5)
Oftentimes we want to re-order the tree according to a phylogenetic rooting.
The phytools package provides a midpoint.root()ing function for re-ordering the tip labels so that the root is at the two longest branchest.
tree_raxml_midpointroot <- phytools::midpoint.root(tree_raxml)
We can compare the order of the original tree vs the midpoint by accessing data stored in the tree object
compare_tips <- data.frame('original'=tree_raxml$tip.label,
'midpoint'=tree_raxml_midpointroot$tip.label)
compare_tips
## original midpoint
## 1 sp|Q9HZP8|FABV_PSEAE sp|O24990|FABI_HELPY
## 2 sp|P24182|ACCC_ECOLI sp|P54616|FABI_BACSU
## 3 sp|O24990|FABI_HELPY sp|Q81GI3|FABI_BACCR
## 4 sp|P54616|FABI_BACSU sp|Q6GI75|FABI_STAAR
## 5 sp|Q81GI3|FABI_BACCR sp|P0AEK4|FABI_ECOLI
## 6 sp|Q6GI75|FABI_STAAR tr|A0A0M3G239|A0A0M3G239_HAEHA
## 7 sp|P0AEK4|FABI_ECOLI tr|A0A1R0ED53|A0A1R0ED53_HAEPA
## 8 tr|A0A0M3G239|A0A0M3G239_HAEHA tr|A0A4Y9GUQ1|A0A4Y9GUQ1_9PAST
## 9 tr|A0A1R0ED53|A0A1R0ED53_HAEPA tr|E6KYQ2|E6KYQ2_9PAST
## 10 tr|A0A4Y9GUQ1|A0A4Y9GUQ1_9PAST tr|B0BP34|B0BP34_ACTPJ
## 11 tr|E6KYQ2|E6KYQ2_9PAST sp|I6Y778|FABG4_MYCTU
## 12 tr|B0BP34|B0BP34_ACTPJ sp|P0AEK2|FABG_ECOLI
## 13 sp|I6Y778|FABG4_MYCTU sp|P71534|MABA_MYCS2
## 14 sp|P0AEK2|FABG_ECOLI sp|P9WGT3|MABA_MYCTU
## 15 sp|P71534|MABA_MYCS2 sp|Q9HZP8|FABV_PSEAE
## 16 sp|P9WGT3|MABA_MYCTU sp|P24182|ACCC_ECOLI
## 17 sp|Q73Q47|FABV_TREDE sp|Q73Q47|FABV_TREDE
## 18 sp|Q2P9J6|FABV_XANOM sp|Q2P9J6|FABV_XANOM
## 19 sp|Q9KRA3|FABV1_VIBCH sp|Q9KRA3|FABV1_VIBCH
## 20 sp|Q62L02|FABV_BURMA sp|Q62L02|FABV_BURMA
We can also compare them by plotting both trees side by side
# tree 1
tree_original_visual <- ggtree(tree_raxml)+
geom_tiplab()+
geom_nodelab(aes(label=label))+
geom_treescale(x=10)
tree_midpoint_visual <- ggtree(tree_raxml_midpointroot)+
geom_tiplab()+
geom_nodelab(aes(label=label))+
geom_treescale(x=10)
ggpubr::ggarrange(tree_original_visual, tree_midpoint_visual, nrow=1)
The tree on the right is far easier to understand. The different gene types almost cluster together in monophyletic branches.
Another option to re-organize the tree is to root the tree to an ancestral gene or to an outlier.
We will try this by using ape::root() on the gene sp|P24182|ACCC_ECOLI and tr|E6KYQ2|E6KYQ2_9PAST
# first outgroup tree
tree_raxml_outgrouproot <- ape::root(tree_raxml ,outgroup='sp|P24182|ACCC_ECOLI')
tree_outgroup_visual <- ggtree(tree_raxml_outgrouproot)+
geom_tiplab()+
geom_nodelab(aes(label=label))+
geom_treescale(x=10)
# second outgroup tree
tree_raxml_outgrouproot2 <- ape::root(tree_raxml ,outgroup='tr|E6KYQ2|E6KYQ2_9PAST')
tree_outgroup2_visual <- ggtree(tree_raxml_outgrouproot2)+
geom_tiplab()+
geom_nodelab(aes(label=label))+
geom_treescale(x=10)
# plot together
ggpubr::ggarrange(tree_outgroup_visual,tree_outgroup2_visual,nrow=1)
The second tree provides a similar ordering that accentuates the divergence of fabV.
Use ape::read.tree() to read in the tree in 'RAxML_bipartitions.raw_enoyl_seqs' and store it to the variable my_tree_raxml.
Use ape::root() to root the my_tree_raxml to the outgroup = 'sp|P0AEK4|FABI_ECOLI' and store it in the variable my_rooted_tree
Visualize the my_tree_raxml using ggtree() with with annotation layers geom_treescale(x=10), geom_tiplab(), and geom_nodelab(aes(label=label))
The ggtree package offers a lot of customization options for your tree plots similar to ggplot2.
We will now go over a variety of them below.
First we will make a tree for a base comparison
tree_raxml <- ape::read.tree('RAxML_bipartitions.raw_enoyl_seqs')
tree_raxml_outgrouproot_compare <- ape::root(tree_raxml ,outgroup='tr|E6KYQ2|E6KYQ2_9PAST')
tree_raxml_outgrouproot_visual <- ggtree(tree_raxml_outgrouproot_compare)+
geom_tiplab()+
geom_nodelab(aes(label=label))+
geom_treescale(x=5)
Making tip labels smaller, making bootstrap values smaller and moving them a bit so they over lap less
tree_compare <- ggtree(tree_raxml_outgrouproot_compare)+
geom_tiplab(size=2)+
geom_nodelab(aes(label=label),nudge_x=-.18,nudge_y=.4,size=2)+
geom_treescale(x=5)
ggpubr::ggarrange(tree_raxml_outgrouproot_visual,tree_compare,nrow=1)
Adding tip points and offsetting the tip labels to accomodate them
tree_compare <- ggtree(tree_raxml_outgrouproot_compare)+
geom_tiplab(size=2,offset=0.1)+
geom_tippoint()+
geom_nodelab(aes(label=label),nudge_x=-.18,nudge_y=.4,size=2)+
geom_treescale(x=7)
ggpubr::ggarrange(tree_raxml_outgrouproot_visual,tree_compare,nrow=1)
Convert the tree into a cladogram
tree_compare <- ggtree(tree_raxml_outgrouproot_compare, branch.length = 'none')+
geom_tiplab(size=2,offset=0.3)+
geom_tippoint()+
geom_nodelab(aes(label=label),nudge_x=-.6,nudge_y=.4,size=2)+
geom_treescale(x=25)
ggpubr::ggarrange(tree_raxml_outgrouproot_visual,tree_compare,nrow=1)
Have all tip labels be right-justified
tree_compare <- ggtree(tree_raxml_outgrouproot_compare)+
geom_tiplab(size=2,align = TRUE,offset=.3)+
geom_tippoint()+
geom_nodelab(aes(label=label),nudge_x=-.18,nudge_y=.4,size=2)+
geom_treescale(x=7)
## Warning: Duplicated aesthetics after name standardisation: size
ggpubr::ggarrange(tree_raxml_outgrouproot_visual,tree_compare,nrow=1)
Cladogram, circular tree
tree_compare <- ggtree(tree_raxml_outgrouproot_compare, layout='circular',branch.length = 'none')+
geom_tiplab(size=2,offset=.01)+
geom_tippoint()+
geom_nodelab(aes(label=label),nudge_x=-.18,nudge_y=.4,size=2)+
geom_treescale(x=20)
ggpubr::ggarrange(tree_raxml_outgrouproot_visual,tree_compare,nrow=1)
Circular tree with branch lengths, all tips aligned to the same branch length
tree_compare <- ggtree(tree_raxml_outgrouproot_compare, layout='circular')+
geom_tiplab(size=2,align=TRUE,offset=.01)+
geom_tippoint()+
geom_nodelab(aes(label=label),nudge_x=-.18,nudge_y=.4,size=2)+
geom_treescale(x=4)
## Warning: Duplicated aesthetics after name standardisation: size
## Warning: Duplicated aesthetics after name standardisation: size
ggpubr::ggarrange(tree_raxml_outgrouproot_visual,tree_compare,nrow=1)
Read in two trees: 'raw_enoyl_seqs.nwk' and 'RAxML_bipartitions.raw_enoyl_seqs' with ape::read.tree(). Store in the variables tree1_fasttree and tree2_raxml
For both tree1_fasttree and tree2_raxml, use ape::root() to root the tree to the leaf using outgroup='tr|E6KYQ2|E6KYQ2_9PAST'
Run ggtree() on both both tree1_fasttree and tree2_raxml and store them in variables tree1_fasttree_visual and tree2_raxml_visual
Use ggpubr::ggarrange() to plot both tree1_fasttree and tree2_raxml side by side
Add the annotation layer scale_x_reverse() to flip the orientation of the tree2_raxml_visual tree
Use ggpubr::ggarrange() to plot the flipped tree2_raxml_visual next to tree1_fasttree_visual
One last thing we’ll touch on are zooming in on specific areas of the tree and collapsing branches.
The nodes of each tree can be visualized
tree_compare <- ggtree(tree_raxml_outgrouproot_compare)+
geom_tiplab(size=2,align=TRUE,offset=.01)+
geom_nodepoint(color='blue',shape=9)+
geom_tippoint()+
geom_nodelab(aes(label=label),nudge_x=-.08,nudge_y=.4,size=2)+
geom_treescale(x=4)
## Warning: Duplicated aesthetics after name standardisation: size
tree_compare
In each ggtree object there exist a list of node numbers. These numbers are different from the bootstrap numbers assigned to a node.
We can list ALL node numbers (including those that belong to leaves) on the tree using geom_text2()
Importantly, we display the values by using aes(label=node)
ggtree(tree_raxml_outgrouproot_compare)+
geom_text2(aes(label=node), hjust=-.3,size=3)
This gives us all node numbers. But we really only want those that are at actual nodes at each branching and not those that are on the leaves. To do so, we use geom_nodelab().
Again, we show node values by using aes(label=node)
ggtree(tree_raxml_outgrouproot_compare)+
geom_nodelab(aes(label=node),hjust=-.3,size=3,color='blue')+
geom_treescale(x=2)
Knowing these node values is useful because now we can use this information to subset the tree and zoom in on a specific portion of it.
tree_raxml_outgrouproot_visual <- ggtree(tree_raxml_outgrouproot_compare)+
geom_nodelab(aes(label=node),hjust=-.3,size=3,color='blue')+
geom_tiplab(size=3)+
geom_treescale(x=7)
zoomin_visual <- viewClade(tree_raxml_outgrouproot_visual, node=25)
ggpubr::ggarrange(tree_raxml_outgrouproot_visual,zoomin_visual)
If we know the structure of the tree, and its nodes, we can clean up the code so that we show bootstraps instead of internal node numbers.
tree_raxml_outgrouproot_visual <- ggtree(tree_raxml_outgrouproot_compare)+
geom_nodelab(aes(label=label),nudge_x=-.15,nudge_y=.4,size=2)+
geom_tiplab(size=3)+
geom_treescale(x=7)
zoomin_visual <- viewClade(tree_raxml_outgrouproot_visual, node=25)
ggpubr::ggarrange(tree_raxml_outgrouproot_visual,zoomin_visual)
Read in the 'RAxML_bipartitions.raw_enoyl_seqs' tree using ape::read.tree(). Store in the variable tree1_raxml. Root the tree to 'tr|E6KYQ2|E6KYQ2_9PAST' using ape::root. Make a tree, tree1_raxml_visual with ggtree() showing each internal label=node using geom_nodelab()
Zoom in on node 35 with viewClade() and store the zoomed in visual in the variable zoomin_visual1
Branches can be collapse()’d at specified nodes. The node can be marked with geom_point2() to more clearly indicate the branch was collapsed.
# collapse, no marker
tree_raxml_outgrouproot_visual2 <- collapse(tree_raxml_outgrouproot_visual, node=36)
tree_raxml_outgrouproot_visual2
# collapse, with marker
tree_raxml_outgrouproot_visual2 <- tree_raxml_outgrouproot_visual2+
geom_point2(aes(subset=(node==36)))
tree_raxml_outgrouproot_visual2
# Multiple nodes collapsed
tree_raxml_outgrouproot_visual3 <- collapse(tree_raxml_outgrouproot_visual2, node=24)+
geom_point2(aes(subset=(node==24)))
tree_raxml_outgrouproot_visual3
Nodes can be collapsed but still show the extent of the phylogenetic diversity that they capture
collapse(tree_raxml_outgrouproot_visual, node=36, 'max')
collapse(tree_raxml_outgrouproot_visual, node=36, 'min')
collapse(tree_raxml_outgrouproot_visual, node=36, 'mixed')
Branches can be highlighted using geom_highight()
tree_raxml_outgrouproot_visual4 <- tree_raxml_outgrouproot_visual+
geom_highlight(35, 'lightblue')
tree_raxml_outgrouproot_visual4
Branches and clades can be annotated using geom_cladelabel()
tree_raxml_outgrouproot_visual+
geom_cladelabel(node=35, label='fabG', align=TRUE,offset = 2, angle=90, vjust=2, hjust=0.5)
## Warning: Ignoring unknown parameters: vjust
One last thing we might want to do is not display any bootstrap value less than 50. This requires altering the tree itself, not the tree visualization object.
tree_raxml <- ape::read.tree('RAxML_bipartitions.raw_enoyl_seqs')
tree_raxml_outgrouproot <- ape::root(tree_raxml ,outgroup='tr|E6KYQ2|E6KYQ2_9PAST')
tree_raxml_outgrouproot$node.label[as.numeric(tree_raxml_outgrouproot$node.label)<50] <- '' # make bootstrap values less than 50 an empty character
## Warning in
## tree_raxml_outgrouproot$node.label[as.numeric(tree_raxml_outgrouproot$node.label)
## < : NAs introduced by coercion
tree_raxml_outgrouproot$node.label[as.character(tree_raxml_outgrouproot$node.label)=='Root'] <- '' #make the word `Root` an empty character
Now we can combine all of the above steps to make an informative plot.
We know our node numbers from previous visualizations but we want to show bootstrap values. We also want bootstrap values less than 50 to not be displayed since we deem them uninformative and cluttering the visualization.
We also want to highlight that fabG has some variation but that one fabG, which is novel, is on a different clade, closer to FabV. We are not too interested in fabV other than it’s location, so we want to collapse it.
Putting it all together.
tree_raxml_outgrouproot_visual6 <- ggtree(tree_raxml_outgrouproot)+
geom_tiplab(size=2,offset=0.1)+
geom_tippoint()+
geom_nodelab(aes(label=label),nudge_x=-.1,nudge_y=.4,size=2)+
geom_treescale(x=7, color=NA)
tree_raxml_outgrouproot_visual6 <- collapse(tree_raxml_outgrouproot_visual6, node=24,'max')
tree_raxml_outgrouproot_visual6 <- tree_raxml_outgrouproot_visual6+
geom_highlight(35, 'lightblue')+
geom_cladelabel(node=35, label='fabG', align=TRUE,offset = 1.2, angle=90, vjust=2, hjust=0.5)+
geom_cladelabel(node=24, label='fabV', align=TRUE,offset = -.2, angle=90, vjust=2, hjust=0.5)+
geom_cladelabel(node=2, label='novel\nfabG', align=TRUE,offset = 1.1, angle=90, vjust=2, hjust=0.5, color='blue')
## Warning: Ignoring unknown parameters: vjust
## Warning: Ignoring unknown parameters: vjust
## Warning: Ignoring unknown parameters: vjust
tree_raxml_outgrouproot_visual6
## Warning: Removed 5 rows containing missing values (geom_point_g_gtree).
There are many more options for manipulating trees available, view them here!
Thanks for coming to Day 1!